Application of elastostatic Green function tensor technique to electrostriction in 

cubic, hexagonal and orthorhombic crystals 



J. Hlinka 1 and E. Klotins 2 

1 Institute of Physics ASCR, Praha, Czech Republic 
1 Institute of Solid State Physics, University of Latvia, Riga, Latvia 
g ! (Dated: February 1, 2008) 

The elastostatic Green function tensor approach, which was recently used to treat electrostric- 
£N| ' tion in numerical simulation of domain structure formation in cubic ferroelectrics, is reviewed and 

extended to the crystals of hexagonal and orthorhombic symmetry. The tensorial kernels appearing 
in the expressions for effective nonlocal interaction of electrostrictive origin are derived explicitly 
and their physical meaning is illustrated on simple examples. It is argued that the bilinear cou- 
pling between the polarization gradients and elastic strain should be systematically included in the 
Ginzburg-Landau free energy expansion of electrostrictive materials. 
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I. INTRODUCTION 



Technological advancements in ferroelectric materials have triggered interest in the kinetics of domain pattern 
formation and its description by the time dependent Ginzburg-Landau model. Unlike the standard thermodynamic 
theory, valid for homogeneous monodomain crystals, the time dependent Ginzburg-Landau model comprises extra 
terms accounting for the contribution of long-range interactions in the free energy and providing a potential to 
simulate kinetics of the domain patterns, macroscopic ferroelectric response, the impact of defects, and gives a new 
insight in the piezoelectric effect. Crucial for this technique, based on a variational derivation of the free-energy 
density with respect to the polarization, is dealing with the elastic field controlled by the inhomogeneous polarization 
due the electrostrictive coupling. Indeed, the static and kinetic properties of the domain structure are substantially 
influenced by elastic strain fields associated with polarization inhomogeneities. For example, it was shown that the 
electrostrictive interactions are critical to the formation of twin domain structure in prototype ferroelectric material 
BaTi0 3 .El 

The standard approach to the problem consists in elimination of the elastic degrees of freedom with the help of 
the mechanical equilibrium conditions. This leads to an effective interaction term _Fhot depending explicitly on the 
order parameter (polarization P(x)) only. Such an effective term is then added to the Landau-Devonshire free energy 
functional instead of the elastic and the electrostrictive terms. A technical drawback of this approach is that the 
resulting effective energy term is nonlinear and nonlocal. Moreover, if the real space integration is preserved, the 
\Q • elastic properties of the medium come in the formula through the anisotropic elastostatic Green functions for which 
' only complicated integral expressions are known. 

Nevertheless, at least for some purposes, the explicit expressions for anisotropic elastostatic Green functions can be 
avoided by expressing the interaction _Fhct in terms of Fourier components of polarization. The Fourier representation 
is particularly convenient in the case of modulated ferroelectrics, where the polarization has a form of a single plane 
waveB, but it is extremely useful even in the case of 3D domain structuresEJEI. In the Fourier representation, the elastic 
properties appear in the expression for F^ c t through a tensorial kernel Bijki(n). This tensorial kernel is 4-th order 
■ tensor angular function comprising all necessary information about the electrostrictive and the elastic properties of 
the material. The elastic properties of the medium are introduced in By/-{(n) solely via the so called elastic Green 
function tensorta, which is a much more simple object than the (real space) elastic Green function itself. 
J-^ ' In some cases, the, exact form of the tensorial kernel B^ki (n) could be reasonably approximated by that of elastically 
^ , isotropic jpediumD'El. Some time ago, however, the components of B.^ki (n) in general cubic crystals-pSwere derived 
k>( ' explicitlyB, and the fully anisotropic i*h e t was then successfully used in realistic 2D and 3D simulationsotm of domain 
rS | structure coarsening in perovskite ferroelectrics. Objective of this work is to the extension of this technique, which 
i may be called " elastostatic Green function tensor technique" , to the crystals of lower symmetry. 

For the sake of clarity, we have introduced the notation and the approach leading to the expressions for effective 
energy contribution F\ lct in section II. The explicit expressions known for cubic crystalsu are generalized to the case of 
orthorhombic and hexagonal symmetries in section III. The section IV is devoted to the basic electrostriction (without 
gradient terms). The .Fhet is expressed in terms of polarization autocorrelation tensor, and the physical meaning of 
Aijki tensor introduced in Ref. |^ is discussed in detail. Finally, the role of gradient terms in systematic expansion of 
electrostrictive energy is elucidated in the Section V. 
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II. ELIMINATION OF ELASTIC DEGREES OF FREEDOM 



The excess Gibbs free energy functional describing an elastically linear ferroelectric in a general polarization and 
stress state can be expressed as a sum of three terms 

F = F {Pi,P itj } + F x {P h P^, Uij } + F 2 { Uij }, (1) 

where Pi, Pij stands for the i— th cartesian component of the polarization field and for its j — th spatial derivative, 
and Uij is the ij— component of the (infinitesimal) strain field. 

The first part Fo{Pi, Pi.j} may be further divided into an integral of the basic local Landau free energy density /l, 
Ginzburg (gradient) energy density /g depending on spatial derivatives of P(x) and the contribution of dipole-dipole 
interaction F^ p 



F {Pi, P i:j } = F dip + J(h+ /g) dx . 



(2) 



The electrostrictive energy Fi{Pi, Pi t j,Uij} can be expressed as an integral over electrostriction density / es , which is 
by definition linear in the elastic strain field Ujj (x) 

Fl{Pi, PiJ : Uij} — j f cs dx, / cs = QijUij, (3) 

where the leading term in expansion of coefficient g%$ 

u. : = .'/,,■: /'•/'...; = q i3 kiPkPi + ••• (4) 

is just given by the usual electrostriction tensor qijki- Note that Einstein summation rule is assumed in the paper. 
The last term, the proper elastic energy, is merely a quadratic function of the elastic strain field 

F 2 {Uij} = J (jPijUij + /cla)dx, /da = ^CijklUijUkl ■ (5) 

In this case the total stress (Jij (x) can be divided in three contributions — thermal stress , describing for example 
the common thermal dilatation, purely elastic stress CijkiUki and the proper electrostrictive tensile stress field 
originating from coupling to the polarization field: 

df 

(Tij(x) = = Pij + CijkiUki - Qiy (6) 

OUij 

In static problems or in dealing with slow processes like domain structure formation, the inhomogeneous elastic strain 
can often be eliminated by means of static equilibrium conditions. The local stress equilibrium condition cr^-j = 
can be considered as a second order partial differential equation for displacement field u(x): 

d 2 u k _ dgij 

^ijklT: — 7^ — — t; — ■ I'J 

OXjOXl OXj 

This condition defines u(r) up to a linear form (homogeneous strain). Let us first consider a macroscopically clamped 
crystal with a large volume V, where the homogeneous strain is zero. In principle, the solution to the Eq. (0) satisfying 

Uij = (Hy(x)) = j Ujj(x)dx = 0, (8) 

can be found using the corresponding anisotropic elastostatic Green function G(x) defined by 

d 2 Gk 

Cijkl Q x _Q^ i = S (x) S im , (9) 

where <5(x) and 8{j are Dirac and Kronecker deltas. 

At the same time, the formal solution for fc =/= Fourier components follows immediately from Eq. (^): 

, / n / , n> ifi(n) • o(k) • k , , 

u(k) = (u(x) exp(-ikx)) = V '-^ ' , (10) 
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where n is the unit vector such that k = fen, 

ffii(k) = (fify(x)exp(-ikx)) (11) 

are the Fourier components of the electcostrictive tensile stress field gij(x) and the Green function tensorB 51(n) is 
the inverse of the Christofell (acousticalQ) tensor T(n): 

(^(n) -1 )^ = T lj (n) = Cijklfikfii- (12) 

Assuming a Born-Karman-like boundary conditions on the volume V , the inverse Fourier transform provides the 
heterogeneous strain field asB 

Uij(x) = ikiUj(k) exp(ikx). (13) 

The k- vectors involved in the summation form a discrete set spread homogeneously over the whole Brillouin zone with 
density V/(2ir) 3 . In fact, only the long wave contributions should be essential (theory assumes smooth inhomogeneity 
or thick enough domain walls) since otherwise the long wavelength elasticity considered in Eq. (|^) is not adequate. By 
inserting the formal solution Eq. (|l^) back to Eq. (Gfenwe obtain the searched effective interaction term F^et = F 1 +F 2 
for a macroscopically clamped system in the formffiffl 



fc^O ijkl 

= -|E fi '5(k)-^(n)-5(-k)-n, (14) 

fc#0 

which does not depend on the elastic strain field any more. The integrand (summand) of Eq. ( |l4| ) is a bilinear form, 
in Fourier transformed tensor components gij (k) . It is thus possible to use the Voigt abbreviated subscript notations 
for Cijki and for g i3 

.9i = .9n, 52 = 522, 53 = .933, 

54 = 523 , 55 = 513 , 56 = 512, (15) 

and rewrite Eq. (jl4|) in a compact form! (hereafter the Greek indices go always from 1 to 6) 

F het = -^^ 5a (k)£^(n) 5/3 (-k). (16) 

The nonlocal character of this effective interaction is better apparent after returning back to the real space : 

Fhct = ~\ Iv ly ~ 9ij (X ° ^dl^dx^ ~ 9lm[x " ] dx ' dx,/ (17) 

where g^ = g^ — (gij) is the heterogenous part of electrostritive field g^. This expression shows that the B a p 
tensorial kernel is actually a Fourier-transformed Hessian of the elastostatic Green function G(x) defined in Eq. (|J) . 
Finally, the total strain field under general macroscopic equilibrium conditions reads 

Mij(x) = Uij + Uij(-x), (18) 

where u,ij(x) is given by Eq. (|l|) and is the homogeneous component defined by the left hand side of Eq. (||) . For 
example, the free sample condition ct^ = leads to the equilibrium value ( see Eq. (||)) 

Utj = Sijkl(Skl - Pki), (19) 

where Sijki = (C~ 1 )ijki is the matrix of elastic compliances. Substitution of Eq. ( |l8| ) back in the original potential in 
Eq. (|l|) provides Fx + F 2 = F hct + F hom where 

V _ 

Fhom = --^{ga-Pa)S a /3(gp-pp) (20) 

and Fhet is just the same as for the case of clamped crystal (Eqs. (Ill), (Ft 
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III. ELASTOSTATIC GREEN FUNCTION TENSOR FOR CUBIC, HEXAGONAL AND 

ORTHORHOMBIC CRYSTALS 

Provided that the explicit dependence of the electrostrictive tensile stress gtj on polarization field appearing in 



the Eq. is known, the effecive energy teem Fh e t can be calculated from Eq. ([14|), (16) or (|17j). Obviously, in 
some caseaS it is worth to tackle the problemliil of ca 



l fttion of the Green function derivatives appearing in rSo. fll'i 
explicitly, while in other cases it is possible to avoid ittlSuu, and use Eq. (|l4| ) or ([l6|). For example, it was showntrH that 
simulations of domain structure coarsening described by the time dependent Ginzburg-Landau equations (including 
the above effective interaction term -Fhet) can be performed entirely in the Fourier space. 

In this paragraph, we will concentrate on properties of the elastic Green function tensor fiy (n) and the 6x6 matrix 
of the tensorial kernel B a p(fi) appearing in the Eq. (|l6|). In order to avoid numerical inversion of the Christofell tensor 
(Eq. ([[2])) at each wave vector direction n, several authors have derived explicit formulas for the cubic symmetry 
Green function tensor fijj(n). Among them, the approach of the Ref. || is the most suitable here since it allows 
generalization to the case of hexagonal and orthorhombic symmetry. The essential step consists in expressing I\j as 
a sum of diagonal part dj(n) and a tensorial square of a convenient real vector v : 

Tij (n) = dj(fi)6ij + ViVj. (21) 

This is trivial in cubic crystals where obviously^ 

Vi = (C u + C44K , (22) 

d i = C 44 + {C n -C 12 + C 44 )hl (23) 

The decomposition is not so straightforward for the crystals of lower symmetry. Nevertheless, for example in the case 
of orthorhombic elastic medium with 

C23 > — C44 , C13 > — C55 , C12 > — C*66 , (24) 

(which is a very weak assumption since practically all known crystals have all the off-diagonal elements C12, C13 and 
C23 positive), the Christofell tensor is given by Eq. @ with 



Vl = Illy 



(C12 


+ C 66 )(C 13 -f 


C55) 


(C23 + C44) 


(C23 


+ C44)(Ci2 4 


Cm) 


(C13 + C55) 


(C 13 


+ C 55 )(C 23 4 


C44) 


(C\2 + Cqq) 



V3=n 3 A j- — —-— , (25) 



di = Cun? + Ce 6 nl + ^5^3 - u? , 

d2 = Ci4nf + (722^2 + C44«3 — v\ , 

da = C^n 2 , +C 44 n 2 2 +C 33 hl-vl (26) 

Obviously, the above decomposition can be used also for hexagonal crystals; it is sufficient to put C55 = C44, C22 = Cn, 
C23 = C%3 and 2C*66 = C11 — C%2- 

For arbitrary crystal symmetry, once the explicit expressions for Uj and di are known, the Green function tensor 
Qij is obtained directly using the Lemma from the Appendix: 

Mfi) = ¥-l?(i+y:#r 1 - (27) 



dj didj dk 



The tensorial kernel B a p then reads 



3 V 2 

^ d k ' 



B a p{ii) =/3 a p-9 a 9p(l+ V^)- 1 , (28) 
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where 

61 = hivi/di, 2 = n 2 v 2 /d 2 , 3 =n 3 v 3 /d 3 , 

4 = n 2 v 3 /d 3 + n 3 v 2 /d 2 , 

05 = n 1 v 3 /d 3 + n 3 vi/di , 

06 = n 2 vi/d 1 + n 1 v 2 /d 2 , 

modifies the eq. (4.16) used in Ref. || for cubic crystals and the components of the {Pap} tensor 
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as a function of di simply coincide with those given previously for cubic crystals in eq. (4.15) of 



TOO Of 



(29) 



(30) 



Let us note that in the rare cases when some of denominators in Eqs. ^2jj) would become zeroO of negative, the 
method works equally well, it is sufficient to modify these equations in order to express the Christofcll matrix in the 
form assumed in Appendix. 



IV. BASIC ELECTROSTRICTION IN CRYSTAL OF ARBITRARY SYMMETRY CLASS 

In this section, we will assume that the electrostrictive tensile stress gij is a bilinear form of polarization components 

g i:j (x) = q ijk iP k (x)Pi(x), (31) 

where qijki is the usual electrostrictive tensor, symmetric both in the first and second pair of indexes. In all crystal 
symmetry classes, at least some of the components are nonzero. Since it is also the lawest parder term in non- 
piezoelectric materials, jnost of the phenomenological models are limited just to that termEMuBtlJ In this case, it is 
convenient to introduced an autocorrelation tensor (k) 



Yy(k) = (P J (x)P,(x)exp(-ikx)) , 
which is nothing else but convolution of corresponding Fourier components of the polarization field 

Y-„-(k) =VJp i (k')P,-(k-k'). 

k' 



The heterogeneous effective energy term F\ lct then readsl 

V 



/"!„, : : --J2Y a (k)A a p{n)Y (-k). 



(32) 



(33) 



(34) 



where A a s(fi) — q a fjBfj~ f (xi)q~ f 5 now depends on both elastic and electrostrictive material constants. 

Let us now assume that the sample contains a single planar domain wall perpendicular to a fixed direction fio. Then, 
we can keep only P(k) and Y Q (k) with k || rio so that A a/ 3(n) = A a p(ho) can be taken in front of the summation 
symbol in Eq. (|||) 



V 



"hot 



A a (3(n ) 



^Y a (k)Y (k)-Y a (0)^(0) 



(35) 



Let us further assume, for example, a 180° domain wall with P\,P 3 — 0, so that only Y 2 (k) contributes. The Einstein 
summation in the above expression then reduces to a single term 



V 



"hot 



A 22 {n ) [(Pi) (Pf) 5 



(36) 
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where () stands for the spatial average as defined in Eq. (g). 

Since the polarization field comes in Eq. ( |36| ) via spatial mean square deviation of P^ , it is apparent that the nonzero 
contributions to Phct comes only from the region in the vicinity of the domain wall. Our example thus allows to give 
a clear interpretation to the A a f}(iio) tensor function. For a fixed domain wall profile, its angular dependence defines, 
how the electrostrictive reduction of domain wall energy varies with domain wall orientation no, and its tensorial 
components distinguish various type of domain walls according to the associated change in polarization direction. 
Obviously, in a well-coarsened domain pattern, where the polarization inhomogeneities are limited to the domain wall 
regions, the Phet become again effectively local functional, but depending on density, type and orientation of domain 
walls. 

Incommensurate structures with modulated polarization represent another transparent case of ID inhomogeneity 
where Eq. ( |35| ) holds. Actually, it was showncl that the P rc i term is essential for explanation of the dielectric anomalies 
of incommensurate ferroelectric NaN02 . In the ideal case of uniaxial sinusoidal modulation with Pi , P3 = and wave 
vector ko || no, there is only a single nonzero pair of Fourier components of polarization {P2(ko), P2( — ko)} so that 

P hct = -tM 2 , 2 (n ) [Pi(k )Pi(-k )] 2 , (37) 

As noted previouslyi, this expression does not depend explicitly on modulation wave vector, and thus it does not 
vanish in ko — > limit. In fact, this observation is not so surprising in the present context, since it follows from the 
fact that the volume ratio between "domain walls" and "domains" is fixed by sinusoidal profile of the modulation. 
Naturally, if one assumes that with decreasing ko the modulation becomes of a more and more "rectangular" shape, 
the "gap" between the energy of homogeneous and modulated ferroelectric would vanish in ko — > limit. 

Finally, let us note that in the case of ID inhomogeneity with a fixed direction no, the first term on the right hand 
side of Eqs. (^) and (^) can be actually interpreted as a local term, merely renormalizing 4-th order terms in the 
Landau-Devonshire potential Pl- At the same time , the second term, although non-local, depends on polarization 
in the same way as the Phom- 



V. GRADIENT ELECTROSTRICTION 



Since the elastostatic Green function tensor technique described in section II was developed for dealing with in- 
homogeneous polarization configurations, it is natural to include in the free energy expansion terms depending on 
spatial derivations of polarization. In principle, consistent free energy expansion may require such terms not only in 
the expansion of Po, but also in the expansion of Pi. Thus, instead of Eq. (|3l|), one may need to assume (in a crystal 
with centrosymmetric paraelectric phase): 

P P x dPk jl dPk dPm a ™ 

9ij = qijkiPkn + njki-Q^- + Si 3 klmn ~Q^r~Q^. — ^ — ■ 

Let us demonstrate the role of this gradient electrostriction terms in the case of uniaxial (Pi, P3 = 0) ferroelectric 
with orthorhombic paraelectric phase. Due to the choice of easy polarization direction and the obvious symmetry 
constraints, all nonzero terms in <? a (x) expansion up to the second order in P = P2 and Pj — dP/dxj can be easily 
enumerated and conveniently expressed using Voigt notation: 

91 = quP 2 + r 12 P, 2 + si 2i P ■ , 

92 = <722P 2 + r 22 P,2 + S 2 2iP 2 t , 

93 = 932P 2 + r 32 P, 2 + S 3 2jPj , 

54 = 7"4 4 P 3 + S442P3P2 , 

55 = S564P1-P3 , 

96 = ^66-P.l + S 6 62PlP,2- (39) 

The Fourier components g Q (k) then read 

ffi(k) = (<?i2 - s l2i kf)Y(k) - ir 12 k 2 P(k) , 

ffa(k) = (922 - s 22i kf )y(k) - ir 22 fc 2 P(k) , 

93(k) = (532 - s 32i kf)Y(k) - ir 32 k 2 P(k) , 

34 (k) = -ir 44 fc 3 P(k) - S442fc 3 /c2^(k) , 

# s (k) = -S5 64 fcifc 3 r(k) , 

56 (k) = -ir 66 fciP(k) - s^fcifeHk) , (40) 
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where we have used autocorrelation tensor component Y(\t) = Y^fk) defined previously in Eq. (p2[). 

The effective interaction F\ ict can be now evaluated from Eq. (|lq). Let us examine the case of ID inhomogeneity 
where only P(k) with k || fii parallel to the crystallographic axis x\ is nonzero. Then, the only nonzero components 
of the B a p{h\) tensorial kernel (Eq. (|2q)) are 



Sn(fii) = ^- , S BB (fii) = ^- , B 66 (ni) = ^- , (41) 

^55 ^66 



and the F^ct then reduces to sum of two terms 



and 



= -|E ^^)!y(fa)y(-fe) , (42) 



F hct2 =~X) ^m)^(-fci)- (43) 



The first term has the same form as the expression in Eq 
generalized A a p tensor dependent also on the modulus of t 



except for the fact that the s a( 3 7 coupling makes the 
re k vector. Therefore, the s Q ^ 7 tensor terms in Eq. J38| ) 



can be safely neglected (unless the problem under study is drastically sensitive to the inhomogeneity lengthscale, as 
for example at lock-in phase transition in type-II incommensurate systemscl). 
The second term given by Eq. (Ea) can be straightforwardly transformed to 



2 



F -«=-2^/ v (^rJ dx ' (44) 

so that it is apparent that this term renormalize the coefficient of the lowest order gradient term in the "Ginzburg 
part" fa of the free-energy expansion (Eq. (|l|)). The experimental studies of bilinear coupling between soft mode 
and acoustic branches by Brillouin and inelastic neutron scattering techniques show that the ryjw coupling term 



in Eq. (38) may indeed cause an important renormalization of the Ginzburg term. Probably, more pronounced 
effects arc expected in crystals with a small Ginzburg term. It is even believed that in some crystals this gradient 
electrostriction compensate the Ginzburg term completely what leads to the appearance of incommensurate modulated 
structural3 : E£l. Unfortunately, in the general 3D case, the effect of the the ryfe; coupling term does not reduce to a 
simple renormalization of coefficients in the Ginzburg free energy and the full anisotropy of B a p(n) tensorial kernel 
should be taken into account. 



VI. CONCLUSION 



The presented elastostatic Green function tensor technique concerns with the simulation of the ferroelectric domain 
pattern being a cutting edge problem both in theory of phase transitions and technological applications. We have 
found that this technique, applied recently to electrostriction in ferroelectrics with a cubic paraelectric phase can be 
straightforwardly generalized to hexagonal and orthorhombic crystals. The contribution of this approach is most valu- 
able for orthorhombic crystals since they are far from isotropy and the closed formulas for elastostatic Green functions 
are known only for a few very special limit casesfl Unfortunately, the method outlined here is not very convenient 
for tetragonal, trigonal and monoclinic symmetries since the v vector used in decomposition of the Christofell matrix 
would have nontrivial angular dependence. We are not aware of any elegant method for inversion of Christofell matrix 
in such cases. 

The essential effect of the nonlocal, nonlinear and anisotropic effective interaction term Fhet consists in reduction 
of domain wall energies in function of their orientation and the associated change of polarization vector. This 
information is conveniently contained in the tensor A a p{p) introduced in Eq. (|34|). Polar diagrams of the A a ^(n) 
tensorial components may thus be quite instructive for understanding the behavior of a particular system. 

Finally, the bilinear coupling between the polarization gradients and elastic strain should not be overlooked in the 
realistic simulations. The values of the corresponding tensorial coefficients ryjy can be determined for example with 
the help Brillouin and inelastic neutron scattering techniques. 
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APPENDIX 



Let us suppose that a finite regular real symmetrical matrix A can be written sum of a diagonal part D and a real 
multiple of a tensorial (dyadic) square v ® v of a real vector v as 

Aij = (D + Av ® v)ij = djSij + XviVj, (45) 

where A is real. Let Wi = Vi/di. Then, the matrix inverse to A reads 

. i 1 Aw ® w , . 

A = D- 1 - , 46 

1 + Av • w y ' 

provided that the right hand side of Eq. ( |46| ) exists. This can be easily proven by multiplication of expression in 
Eq. © and Eq. @. 

In the case of orthorhombic, hexagonal and cubic crystals, the above result allows to find compact explicit expressions 
for the corresponding elastostatic Green function tensors. For example in Eqs. ( plf ) and (P7|), we have applied the 
Eq. © with A = 1. 
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